 // MyTesTSuite.h
#include <iostream>             
#include <iterator>
#include <cxxtest/TestSuite.h>
#include "smatrix2.hpp"
#include "eliminate_eq.hpp"


#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>



using namespace std;

typedef double coef_t;
using namespace smatrix2;
typedef sparse_matrix2<coef_t> SMatrix;

//using namespace boost::numeric::ublas;
//typedef matrix<coef_t> SMatrix;

typedef vector<coef_t> Vector;

double Round(double v, int decimals){
    int sgn = (int)(v / fabs(v));
    double eps = sgn*pow10(decimals);
    double rmr = fmod( v, eps );
    v -= rmr;
    if( sgn*rmr - sgn*(eps / 2)>-std::numeric_limits<double>::epsilon() )
	v += eps;
    return v;
}	//double Round(double v, int decimals)

class TestF_And:  public CxxTest::TestSuite 
{
	public:

        static TestF_And *createSuite() { return new TestF_And; }
        static void destroySuite( TestF_And *suite ) { delete suite; }
			
	TestF_And(){
	}
	virtual ~TestF_And(){ 
	};
// ++++++++++++++++++++++++++++ Suite ++++++++++++++++++++++++++++++++++++++++
	
	void testAlive( void )
	{
	    TS_TRACE( "Suite started" );
	    TS_ASSERT_EQUALS(2+2,4);
	    cout<<"\nepsilon is: "<<std::numeric_limits<double>::epsilon();
	    //TS_ASSERT_EQUALS(Round(2.3454, -1), 2.3); 
	    //TS_ASSERT_EQUALS(Round(2.3454, -3), 2.345); 
	    
	}
	/// x + y  = 2
	/// x - y  = 0
	
	void test_01(){
	    SMatrix sm(2,2);
	    coef_t l1[]= {1, 1};
	    coef_t l2[]= {1,-1};
	    coef_t c0[]= {2, 0};
	    std::copy(l1, l1+sizeof(l1)/sizeof(coef_t), sm.begind2(0));
	    std::copy(l2, l2+2, sm.begind2(1));
	    Vector c(2);
	    std::copy(c0, c0+2, c.begin());
	    cout<<"\nInput matrix"<< sm<<endl;
	    cout<<"\nInput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    

	    EqSolve <SMatrix, Vector> solver (sm,c);
	    solver.Solve(0,0);
	    
	    cout<<"\nSolved 0"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    TS_ASSERT_EQUALS((int)Round(c[1]/sm(1, 1), -1 ), 1); // = 1
	}
	
	void test_02(){
	    SMatrix sm(3,3);
	    coef_t l1[]= { 3, -2, 1};
	    coef_t l2[]= {-1,  1,  2};
	    coef_t l3[]= { 1, -1, -1};
	    
	    coef_t c0[]= {21, -1, 4};
	    
	    std::copy(l1, l1+sizeof(l1)/sizeof(coef_t), sm.begind2(0));
	    std::copy(l2, l2+sizeof(l2)/sizeof(coef_t), sm.begind2(1));
	    std::copy(l3, l3+sizeof(l2)/sizeof(coef_t), sm.begind2(2));
	    
	    Vector c(sm.size1());
	    std::copy(c0, c0+sizeof(c0)/sizeof(coef_t), c.begin());
	    
	    cout<<"\nInput matrix"<< sm<<endl;
	    cout<<"\nInput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    SMatrix csm = sm;
	    Vector cc = c;
	    
	    //test1
	    
	    EqSolve <SMatrix, Vector> solver (sm,c);
	    solver.Solve(0,0);
	    
	    cout<<"\nSolved 0"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    
	    solver.Solve(1,1);
	    
	    cout<<"\nSolved 1"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    TS_ASSERT_EQUALS((int)Round(c[2]/sm(2, 2), -1 ), 3);
	    if (1) {
	    sm = csm;
	    c = cc;
	    
	    // test 2
	    cout<<"\n===============Test 2====================";
	    	    
//	    EqSolve <SMatrix, Vector, coef_t> solver (sm,c);
	    solver.Solve(0,0);
	    
	    cout<<"\nSolved 0"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    
	    solver.Solve(2,2);
	    
	    cout<<"\nSolved 1"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    TS_ASSERT_EQUALS((int)Round(c[1]/sm(1, 1), -1 ), -3);
	    }
	    if (1){
	    sm = csm;
	    c = cc;
	    
	    // test 3
	    cout<<"\n===============Test 3====================";
	    
	    cout<<"\nInput matrix"<< sm<<endl;
	    cout<<"\nInput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    	    
//	    EqSolve <SMatrix, Vector, coef_t> solver (sm,c);
	    solver.Solve(1,1);
	    
	    cout<<"\nSolved 0"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    
	    solver.Solve(2,2);
	    
	    cout<<"\nSolved 1"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    TS_ASSERT_EQUALS((int)Round(c[0]/sm(0, 0), -1 ), 4);
	    
	    solver.Solve(0,0);
	    
	    cout<<"\nSolved 2"<< sm<<endl;
	    cout<<"\nOutput vector ";
	    std::copy(c.begin(), c.end(), std::ostream_iterator<coef_t>(cout, " "));
	    
	    }
	}
	
};
